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Abstract 

As two-dimensional fluid shells, lipid bilayer membranes resist bending and stretching but are unable to 
sustain shear stresses. This property gives membranes the ability to adopt dramatic shape changes. In this paper, 
a finite element model is developed to study static equilibrium mechanics of membranes. In particular, a viscous 
regularization method is proposed to stabilize tangential mesh deformations and improve the convergence rate of 
nonlinear solvers. The Augmented Lagrangian method is used to enforce global constraints on area and volume 
during membrane deformations. As a validation of the method, equilibrium shapes for a shape-phase diagram 
of lipid bilayer vesicle are calculated. These numerical techniques are also shown to be useful for simulations 
of three-dimensional large-deformation problems: the formation of tethers (long tube-like exetensions); and 
Ginzburg-Landau phase separation of a two-lipid-component vesicle. To deal with the large mesh distortions of 
the two-phase model, modification of vicous regularization is explored to achieve r-adaptive mesh optimization. 

1 Introduction 

Lipid membranes are a critical part of life because they serve as a barrier to separate the contents of the cell from 
the external world. Lipid molecules are composed of a hydrophilic headgroup and two hydrophobic hydrocarbon 
chains [1], and will form a bilayer structure spontaneously by the hydrophobic effect when introduced into water 
in sufficient concentration. Though the cell membrane has more complex structure, being littered with all kinds 
of proteins that serve as selective receptors, channels, and pumps, in this paper we will focus on closed spherical 
pure lipid bilayer membranes, i.e., vesicles. 

Common experience reveals that it is much easier to bend a thin plate than to stretch it (a good example is 
a sheet of paper). This is also true for lipid bilayer membranes, except that there's no shear force because of 
the fluid property of membranes. The mechanical energy of a lipid bilayer has three major contributors; bend- 
ing (curvature) of each monolayer; area or in-plane expansion and contraction of each monolayer; and osmotic 
pressure. Because the last two energy scales are much larger than the first one (by several orders of magnitude) 
[2], they effectively place constraints on the total surface area and enclosed volume of the bilayer membrane on 
experimental time scales (up to at least one hour). Thus, the mechanically interesting energy arises from bending 
of the membrane. 

Canham [3], Helfrich [4] and Evans [5] pioneered the development of the lowest-order bending energy theory, 
often referred to as the spontaneous curvature model, in which energy is a quadratic function of the principle 
curvatures and the intrinsic or spontaneous curvature of the surface. Incremental improvements to this model 
include the bilayer couple model [6, 7], which imposes the hard constraint on the area difference of the two 
monolayers, and the area-difference-elasticity (ADE) model [8-10], which adds a non-local curvature energy 
term representing an elastic penalty on the area difference. 

The equations of equilibrium for the spontaneous curvature model, first calculated by Jenkins [11, 12], are 
difficult to solve being highly nonlinear fourth-order PDEs. The most common approach to modeling membrane 
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mechanics numerically has been to discretize a vesicle surface by a triangle mesh, and approximate the curva- 
ture along mesh edges with finite-difference (FD) operators. Starting from some suitable initial shape, the FD 
approximation of curvature energy can be summed on the triangulation, and an adjacent local minimum then 
can be found by downhill minimization (often via a conjugate gradient algorithm) [13-16]. Another mesh-based 
approach, the finite element method (FEM), was also recently applied to the study of membrane mechanics by 
Feng and Klug [17], using -conforming triangular subdivision-surfaces elements to approximate the membrane 
curvature energy. 

One feature shared in common among these mesh-based methods is the need for stabilization of mesh vertex 
motions tangential to the discretized surface. This issue arises as a fundamental consequence of the use of a mesh 
for explicit coordinate parameterization of the geometry of a fluid membrane having no physically meaningful ref- 
erence configuration. As pointed out previously [e.g., 18, 19], the dependence of the curvature energy functional 
on the surface position map is invariant upon changes in parameterization. Physically this implies that in-plane 
dilatational and shear modes of mesh deformation carry no energy cost, and therefore, no stiffness. The addition 
of an artificial in-plane stiffness, for example by placing Hookean springs along the edges of a triangular mesh 
[18] does indeed stabilize these motions, but in doing so also changes the physics, yielding a model for a poly- 
merized (rather than fluid) membrane. In order to allow for fluid-like diffusion of membrane vertices, a number 
of researchers have used a dynamic triangulation approach, wherein the edge connecting a pair of adjacent trian- 
gles along one of the diagonals between the four associated vertices is swapped for the other diagonal. Within a 
Monte Carlo simulation framework this method yields mean-square vertex displacements that are consistent with 
microscopic diffusion [20], and this approach has been shown to produce an effective viscosity that increases as 
the edge swapping rate decreases [21]. However, it is unclear whether the dynamic triangulation approach can 
enable unphysical in-plane forces to fully relax to zero in a zero-temperature energy minimization context as is 
adopted here. 

Alternatively, as shown in the finite element (FE) context [17], tangential vertex motions may be suppressed 
partially by enforcing the incompressibility of the membrane as a local (rather than global) area constraint. This 
approach also allows for diffusion of vertices; however, local enforcement of incompressibility can lead to severe 
distortion of elements in the mesh, and therefore hinders the simulation of large vesicle deformations. Further- 
more, local incompressibility does not completely suppress spurious modes, and though these degeneracies do 
not prevent simulation of unforced vesicle equilibrium, they can lead to catastrophic numerical instabilities when 
externally applied forces are introduced. 

Notably, these issues may be avoided entirely through the development of meshless numerical methods for 
membrane mechanics. Examples include, Ritz methods with global basis functions (e.g., spherical harmonics) 
[3, 22, 23], phase-field methods [24, 25], and moving-least squares approximation [26]. Yet, these approaches are 
not without their own limitations (e.g., aliasing, difficulty with application of external forces). 

In this paper we propose a viscous regularization technique to stabilize tangential motions of nodes in a FE 
membrane model while enforcing incompressibility as a global constraint. We demonstrate the computational 
efficiency and effectiveness of this approach by comparing simulation times with and without regularization for 
shape transitions previously computed in [ 1 7] . Secondly, we examine the efficiency gained by enforcing the global 
constraints on membrane area and volume with an augmented Lagrangian approach instead of the previous penalty 
approach. Lastly, we apply the regularization and constraint methods to the simulation of two membrane shape 
change problems involving large deformations and the application of external forces. The first of these problems 
is the simulation of the tether instability in a vesicle under tension between two opposing point forces; the second 
is the simulation of separation and domain formation in a two-lipid-phase vesicle. In the latter we demonstrate 
how the viscous regularization technique can be slightly modified to formulate an r-adaptive remeshing method, 
wherein nodes "flow" on the surface of the vesicle in such a way as to avoid element distortion. 

The outline of this paper is as follows; Section 2 briefly introduces the FEM formulation for bilayer membrane 
mechanics along with artificial viscosity mesh stabilization and augmented Lagrangian constraint enforcement; 
Section 3 shows two applications: tether formation and lipid phase separation, based on the methods described in 
Section 2; Section 4 concludes with discussions of results and future applications. 
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2 Methods 



2.1 Lipid bilayer mechanics and finite element approximation 

We begin with a brief review of the mechanics of bilayer membranes and FEM approximation we use. For further 
details, the reader is refer to the paper [17]. 



Z 




Figure 1 : Geometry of a surface 



Membrane kinematics. The bilayer membrane is described as a two-dimensional surface embedded in 
three-dimensional space (Fig.l), parameterized by curvilinear coordinates {s^, s^}, such that its position is given 
by the map a; : ^ M^. With standard definitions from differential geometry [27, 28] we can span the surface 
tangent plane with both (co variant) basis vectors = = x^a and dual (contravariant) basis vectors a" 
defined such that a" • = ^p- The co variant and contravariant surface metric tensors are then 

o-ap = a,a ■ o,fj, and a"'^ = a" ■ a'^ , (1) 



and the determinant of the covariant metric tensor will be denoted 



a = det Qa/s . (2) 

The normal to the surface is 

ai X 02 ai X 02 
d = a3 = - - = (3) 

\ai X a2\ V« 

The curvature tensor B is defined by its covariant components 

bap = -d a ■ ap = d- X,al3 = d ■ aa^0. (4) 

The mean curvature is one half of the trace of the curvature tensor 

H^\bl = ]^a^Pbai3 - -\a''^{d^c. ■ ap) = -^a" • d „ (5) 

(where a"^ is the contravariant metric tensor, defined such that a^^Uap = S^), and the Gaussian curvature is the 
determinant of the curvature tensor 

K = dct B = dct bap. (6) 



Lipid bilayer mechanics. We describe the energetics of the membrane by the Helfrich model [4], which as- 
sumes a strain energy of the form 

E[x\^ ( ]-K,c{2H -Cof^/^(fs+ ( JCoKV^d^s (7) 
Jm 2 Jm 
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where Kc is the bending modulus and TCg is the Gaussian curvature modulus. By the Gauss-Bonnet theorem 
[27], the integral of Gaussian curvature is a topological constant /^^ K^Jad?s — 47r(l — g), with g being the 
genus, i.e. the number of handles, and thus can be neglected. 

The weak form of equilibrium for the membrane can be obtained in general by the principle of virtual work, 
and for the case of conservative loads by minimization of total potential energy. The later dictates that the total 
potential energy by stationary with respect to any arbitrary admissible surface variation 5x 

5Ii = SE- SW""' = 0. (8) 

Here 5E is the first variation of the membrane bending energy, and SW^' = • 5x^/a is virtual work done 

by conservative external forces f^^^. A straightforward calculation [17] gives the first variation of the total energy 
as 

5n= I [n" • (5a„ + • Sd.o, - • 5x\ y/^ds^ds^, (9) 

J M 

where we have defined stress resultants n" and moment resultants m" as 

n" =/Cc(2i?-Co)a"''d^ + /Cc^(2i?-Co)^a" m" = -A:c(2i/ - Co)a". (10) 

Enforcing constraints: augmented Lagrangian method. Admissibility requirements on trial functions x and 
variations 5x include the satisfaction of any active constraints, such as the aforementioned constraints on total 
surface area and enclosed volume. Here we will enforce these constraints with the augmented Lagrangian (AL) 
approach (see, e.g., [29]). The AL method may be thought of as a hybrid between penalty and Lagrange multiplier 
methods. The basic idea of AL is to solve iteratively for a Lagrange multiplier, computing multiplier updates 
from a penalty term. To enforce constraints on both area and volume of a membrane, we establish a sequence of 
modified energy funcitonals, the nth of these taking the the form /" = 11 + where is a constraint energy 
term 

Here A and V are the specified surface area and enclosed volume of the membrane, and /iy are penalty 
parameters (large and positive), and a" and p" are tension and pressure multiplier estimates for the rtth iteration. 
Minimization of the modified energy (holding multiplier estimates fixed) yields 

= - p"+^5V + a^+^5A = 0, 

where = p" — iiv{Vn — V) and = a" + ^A{An — A) are the updated multiplier estimates. Iter- 

ation of minimization followed by multiplier updates is continued until constraints are satisfied to within some 
preselected tolerance, TOL, as shown below in Algorithm 1 . In this way the modified energy converges to the 
pure Lagrange-multiplier constrained functional, with the added benefit of avoiding the associated saddle-point 
problem, retaining a minimization structure which is convenient for nonlinear optimization algorithms. 

Whereas pure penalty methods require very large penalty parameters for accurate constraint enforcement, the 
AL iterative updates can achieve accuracy with much smaller penalty terms. In practice this is an important ad- 
vantage, since when the penalty parameters jiy and jiA become large, numerical minimization becomes difficult 
as the Hessian V^E (or stiffness matrix) becomes quite ill conditioned near the minimizer. This property makes 
minimization algorithms like quasi-Newton and conjugated gradient perform poorly, as finding the search direc- 
tions becomes difficult [29]. However, small penalty parameters can produce a large number of AL iterations for 
convergence. Hence, it is common in practice to incrementally increase penalty parameters by some factor, FAC, 
after each AL iteration process to achieve faster convergence. These penalty parameter updates are also included 
in Algorithm 1 . 

Finite element approximation. A FE approximation is introduced by replacing the field x with the approxi- 
mated field Xh defined by 

N 

Xh{s\s')=Y,XaN''{s\s^) (11) 

a = l 
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Algorithm 1 Augmented Lagrangian method to enforce global area and volume constraints. 

Initialize: p", a", /iy, /i^, 

Set n=0 

repeat 

Minimize /" = H + - Vf - p^V + ^{A- Af + cJ^ A 

Updates: 

^i!^^ = My X FAC , ^i!^^^ = Ml x FAC 
n <— n + 1 

until |y" - 1/| < TOL and |A" ~ A\ < TOL 



where the Af°(s^, s^), a = 1,...,7V are shape functions of the FE mesh, and their coefficients, Xa are the 
positions of the nodal control vertices. Introducing this approximation into the modified energy funcitonal upon 
minimization leads to a set of discrete approximate equilibrium equations 



ST + ST - ST = 0. 

Here are the internal nodal forces due to bending of the membrane. 



/int 
a 



M 



n" • h m" 

OXa 



dd 

dXa 



(12) 



(13) 



/™" are the constraint nodal forces due to the pressure and tension that are conjugate to the constrained volume 
and area. 



" dXa dXa ' 

and fT are the external nodal forces, due to the appUcation of distributed loads on the surface, 



ST= I S'^'N-'T^Ss. 

JM 



(14) 



(15) 



Note that the integrands of the global expressions for internal and constraint forces are described in more explicit 
detail in [17]. Following that work, we again employ -conforming subdivision surface shape functions [30, 31] 
along with second-order (three-point) Gaussian quadrature for the computation of element integrals. 



2.2 Viscous regularization of tangential mesh deformation 

In the curvature model, the energy is determined by the mean curvature which is a parameterization-independent 
property of the surface shape, and thus is not sensitive to in-plane dilatational or shearing deformations of the 
surface FE mesh. Much like physical lipid molecules, FE nodes can flow freely on the deformed surface. As 
discussed in [17], this fact is manifested in the appearance of degenerate, zero-stiffness, zero-energy modes. Here 
we discuss the implementation of an artificial viscosity method designed to numerically eliminate these degenerate 
modes. 

For solid shells having both reference and deformed configurations, in-plane deformations (dilatation and 
shearing) can thus be expressed locally in terms of first derivatives of the surface position maps of these two 
configurations. In curvature model,a well-defined reference configuration does not exist since the energy is only 
related to the deformed shape. The basic ingredients for stabilization of these tangential modes are the introduc- 
tion of a reference configuration and an energy term elastically penalizing in-plane deformation away from this 
reference state. However, to retain the physics of the original model, the addition of any in-plane elastic energy 
must result in a variational problem possessing the same minimizing solution as the original problem. In other 
words, the artificial in-plane energy must attain a value of zero when the entire model is in equilibrium. To de- 
sign an algorithm that achieves these goals, we define a sequence of variational problems, minimizing a modified 
energy functional 

/" = \i[x] + rTx] + I'Tx-, (16) 
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where the reference configuration X" for the nth iteration is the deformed solution a;""^ of the previous iteration. 
The form of the regularization energy r^s[x; X] can be chosen such that it vanishes when x = X, to ensure that 
solutions a;" converge to minimizers of the original unregularized problem with increasing n. This regularization 
method is outUned below in Algorithm 2. 



Algorithm 2 Viscous regularization via reference updates. 
Set = initial shape. 
Set n = 
repeat 

Minimize /" = n[a;]+/™"[a;]+/'''=s[a;;X"] solution, a; 

Update reference: set X"+^ = a;" 
n <— ?i + 1 
until/"=s[a;";X"] < TOL 



Qualitatively, assignment of the reference configuration for each iteration to be the current configuration of the 
previous iteration results in a type of algorithmic viscosity, producing forces that resist the motion of nodes away 
from their position at each previous iteration. The quantitative details of this viscosity depend on the particular 
form chosen for the in-plane regularization energy, 7"^^. Here we give two example forms, the first derived from 
planar continuum elasticity theory and the second representing the mesh as a network of viscous dashpot elements. 

Continuum elastic regularization energy. Here we treat the in-plane deformation response for each regulariza- 
tion iteration as that of a two-dimensional solid membrane. This local response can be modeled via a hyperelastic 
strain energy density, w{F), which is a function of the surface deformation gradient 

F = a„®A", (17) 

where A" are the dual basis vectors on the reference surface, i.e.. A" ■ Ap = 5'^, where A^ = X'^. Thus the 
regularization energy becomes 

/'"=s[a;;X"] /" w{F)^/A<fs. (18) 

JM 

To preserve objectivity, the strain energy is a function of F through implicit dependence on the invariants of 
the surface-Right-Cauchy-Green deformation tensor C = • F = ao,pA°' ® A^ [32]. As C is a rank-2 tensor, 
the two non-zero principal invariants are 

h = tr(C) = a^'^aap (19a) 
1-2 = \{HC)f - tr(C2)} = ^{{a^Pa^pf - a'^^a^"' a^pa^,} = j\ (19b) 

The strain energy density is thus a function of these two invariants 

As a specific example, consider a strain energy function that decouples the dilatational, and shear responses, as 
used by Evans and Skalak [33] to model the red blood cell cytoskeleton 

k, /tr(C) \ 

area change ^f^^^^. 

Here k and i^i are stretching and shear moduli, respectively. It should be carefully noted that although we follow 
here the formalism of solid mechanics, the reference configuration X" is not permanent as for a soUd; rather the 
reference configuration is iteratively updated so that the resulting stresses may relax to zero. 
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Dashpot regularization energy. The viscous character of our proposed scheme is much more obvious when we 
compose the regularization energy of contributions from Hookean springs placed along all element edges, namely, 

jreg^ ^ (20) 
edge ab 

where iab = \xa — xi,\ and Lab = l-X^" — XJ^\ are the lengths of the edge connecting mesh vertices a and b in the 
deformed and current configurations, respectively. Differentiating this energy, the corresponding force on a node 
a from the spring connecting it along an edge to node b can be obtained as 

fab = k{iab - Lab)nab 

where Uab is the unit vector pointing from node a to node b. Recalling that the reference configuration for the nth 
iteration is the same as the deformed configuration of the n — 1th iteration, the magnitude of this force can also 
be written as 

fab = k{£ab ~ Ch )■ 

This is easily identified as the backward-Euler time-discretization of the force-velocity relation for a viscous 
dashpot 

, d£ab 



Thus, iterative reference updates of the form X" = x"^^ have the effect of converting a network of springs into 
a network of dashpots, clearly revealing the viscous character of the regularization scheme. 

We have numerically implemented both the continuum elastic and dashpot regularization described here, and 
although both forms are effective in practice we have preferred the dashpot approach for its simplicity, efficiency, 
and robustness. The remainder of the paper focuses on the use of this second approach, demonstrating its effec- 
tiveness in application. 



3 Applications 

3.1 Shape vs. reduced volume 

Even in the absence of any externally applied loads, the two constraints on area and volume cause vesicles to 
transition among a variety of interesting equilibrium shapes. Here some of the calculations performed in [17] of 
the equilibrium shapes for different reduced volumes are repeated, as a first demonstration of the effectiveness of 
viscous regularization. 

Reduced volume i/ is a geometrical quantity defined as 

V 

(21) 



(4V3)i?o3' 



where Rq = \J AjAn is the radius of a sphere with the area A of the vesicle. Reduced volume is then written as 



^3/2 



(22) 



The reduced volume is the ratio of the current volume of the vesicle and the maximum volume that the current 
total area of vesicle can ensphere. For a spherical vesicle, the reduced volume ly ~ 1; a vesicle of any other shape 
has < < 1. 

To compute the following results, the spontaneous curvature model is used with Co = 0. The modified 
energy is computed with loop subdivision shell elements and second-order (three-point) Gaussian quadrature, and 
minimized with the quasi-newton L-BFGS-B solver [34—36]. 
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Computational cost with and without the viscous regularization 





Total iterations (initial shape v = 0.9) 


Total iterations (v = 0.9 0.8) 





35,950 


603,858 


0.5 


2,304 


287,825 


1 


1,821 


9,482 


10 


2,122 


12,480 


100 


4,776 


79,038 



Table 1: Viscous regularization improves the convergence rate of L-BFGS-B minimization. The first row 
kE? /JCc = 0, shows the results without the viscous regularization, which are identical to the approach used 
in [17]. (k: spring constant; R: average radius of the vesicle; ICc- bending modulus.) 

Viscous regularization. As a first assessment of the benefit of regularization, results are compared with the 
simulations done in [17], in which local area and global volume constraints were performed by penalty method 
instead of AL method. First, the same calculation of [17] is repeated; then the viscous regularization is added, 
with same kind of constraints (local area and global volume constraint) and penalty parameters (/i^ = lQ*R^/ICc 
for local area constraint and ha = 5 x W^R^/ICc for global volume constraint). 

The calculation starts from an initial ellipsoid shape which has a reduced volume i' = 0.914 (Fig. 2). In 
the calculation, the area is fixed at its initial value and the volume is reduced in order to satisfy the constraint 
on ly. For each simulation, violation of the volume constraint subjects the vesicle to a large pressure according 
to the penalty term in the functional. The energy is then relaxed by L-BFGS-B minimization and result in the 
equilibrium shapes. The iteration of of reference updates in Algorithm 2 is continued until the regularization 
energy is sufficiently small, F'^^/I < 10^^. In all the simulations, the same mesh, made up of 642 vertex nodes 
and 1280 elements, is used. 

The resulting equilibrium shapes for v = 0.9 and ly = 0.8 are shown in Fig. 2. Starting from the initial 
shape, the equilibrium shape for v = 0.9 is computed by minimizing the energy; then from the resulting v = 0.9 
shape, setting u ~ 0.8, the equilibrium shape for u = 0.8 is computed. The computational cost with and 
without the viscous regularization is listed in Table 1. As can be seen, the convergence rate is highly improved 
(almost two orders of magnitude faster) with the viscous regularization while the resulting shapes are equivalent. 
For different choices of spring constant k, the computational cost also varies. The computational cost has two 
contributions: one is the iteration number for each minimization; the other is the number of reference updates 
required to satisfy the convergence criterion F'^" /I < 10"^. These both depend on k. For each minimization, the 
larger k is, the smaller the iteration number will be. While for the number of reference updates, it is opposite: 
the larger k is, the more reference updates needed. For example, to get the equilibrium shape ly ~ 0.9 from the 
initial shape, fci?^//Cc=l requires 2 reference updates, each of which costs « 1000 iterations for minimization; 
while for kR^ /ICc=^00, there are 20 reference updates each costing « 250 minimization iterations. In this case 
kR^ / K.c=^ works the best, but the optimal k depends on the specific problem. In the later sections on tether 
formation, a much larger k (/ci?^//Cc'=1000) is used. 

Augmented Lagrangian constraint enforcement. From the results described above, viscous regularization is 
shown to be able to heavily lower the computational cost when a penalty method is used to enforce the constraints 
on area and volume. However, regularization also eliminates the need for local enforcement of incompressibility. 
Global constraints on area and volume can be easily implemented via the augmented Lagrangian (AL) method 
which is more efficient than the previous penalty method. Here, the shape change from the initial ellipsoid shape 
to the equilibrium shape of v ~ 0.9 is used to compare the penalty method with the AL method. In this test 
global area and global volume constraints are carried out first by the penalty method with a range of penalty 
parameters, and secondly with the AL method. Viscous regularization is used for both the penalty method and 
the AL method. Iteration of reference updates is continued until the regularization energy is sufficiently small, 
F'^^/I < 2.0 X 10"^ 

The regularization spring constant k is set to be k = IQJCc/R^, where R is the average radius of the vesicle 
and ICc is the bending modulus. For the AL method, the penalty parameters are initialized to be a fairly small 
number (hv = 10*, ha = 10**), and are then increased by a factor of 2 for each of the following minimizations. 
The minimization continues until the constraints on area and volume are satisfied to within a tolerance and the 
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(a) Initial shape, v = 0.914 (not in equilibrium) 




(b) u = 0.9 without regularization. (c) u = 0.9 with regularization. 




(d) u = 0.8 without regularization. 



(e) u = 0.9 with regularization. 



Figure 2: Limit surfaces and control meshes of equilibrium shapes for u = 0.9 and v = 0.8. 



regularization energy is sufficiently small. Viscous regularization reference updates are included with AL mul- 
tiplier updates in a single iteration loop. This hybrid regularization-AL algorithm, shown in Algorithm 3, is a 
combination of separate Algorithms 1 and 2. 

As the Table 2 shows, the AL method reduces the computational cost significantly. This is especially true 
when high accuracy of the constraints is desired, in which case the penalty method requires extremely large 
parameters, which lead to conditioning problems that impede convergence of the nonlinear solver. Indeed, for 
penalty parameter > 10^, L-BFGS-B iterations diverge. In contrast, for the AL method to achieve high accuracy 
penalty parameters need not be very large [29] . 

3.2 Tether formation 

A point-force acting on lipid membranes can pull out a long narrow tube commonly called a tether. This can 
be done by using micropipettes [e.g., 23], optical tweezers [e.g., 37], or even growing microtubules inside the 
vesicle [37]. The mechanical reason for formation of tethers Ues in the lack of shearing modulus for membranes. 
Elongating in one direction and contracting in the other to such a spectacular way like tethers mechanically means 
extremely large shear deformations [38^1]. 

Since tether simulation involves very large deformations, the triangles in the finite element mesh are subject to 
severe distortions. In practice, as elements become more distorted, the zero-energy tangential modes can actually 
become numerically unstable (Fig. 3). Viscous regularization has to be added in order to suppress these zero- 
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Algorithm 3 Hybrid algorithm combining AL constraint enforcement with viscous regularization. 

Initialize: p", a", /iy, /i^, = initial shape, v = 0.9 (the specified reduced volume). 

Set n = 

repeat 

Minimize /" = n[a;] + ^^(1^ - - + ^(A - A)^ + + /'""^[a;; X"] 
— > solution, a;" 
Update reference: set = a;" 

if — > TOLi then 
AL Updates: 

= - My (K - V) , = a" + MK^n - ^) 

^^+1 = yj^ X FAC , li!^^^ = Ma X FAC 
end if 
n <— 71 + 1 

until |z/„ - i?|/z^ < TOLi and r'=g[a;"; X"] < TOL2 



Penalty method vs. AL method 



Accuracy \v - 0.9|/0.9 


10"^ 


10"^ 


10""* 


10"^ 


10"'' 


10-' 


io-« 


10-*^ 


10-^*^ 


Penalty parameter 


10^ 


10^ 


10*^ 


10'^ 


10« 


n/a 


n/a 


n/a 


n/a 


Iterations (penalty method) 


379 


436 


1032 


3455 


9231 


n/a 


n/a 


n/a 


n/a 


Iterations (AL method) 


366 


425 


500 


515 


616 


743 


916 


981 


1201 



Table 2: Computational cost of the penalty method and the AL method. To achieve the same accuracy, the AL 
method requires fewer total iterations compared to the penalty method. More importantly, for extremely high 
accuracy (< 10~^), the L-BFGS-B minimization diverges with the penalty method, while the AL method still 
converges. 

energy modes. Furthermore, the critical force to pull out a tether is very sensitive to pressure and surface tension. 
Numerically, this necessitates highly accurate enforcement of the volume and area constraints. For a penalty 
method this implies very large penalty parameters, which lead to conditioning problems (e.g.. Table 2). For this 
reason, here the augmented Lagrangian method is applied. 




Figure 3: The equilibrium shape of a vesicle of reduced volume v = 0.9 with external forces (sa 1 pN) applied 
at the two ends (forces not shown), mesh without (left figure) and with (right figure) the viscous regularization. 
Note that the unstabilized mesh is subject to element distortion even at small applied load. 

Starting from an initial equilibrium shape (prolate), tether development is simulated by incrementally displac- 
ing nodes at the tips of the vesicle, and performing energy minimization resulting in the equilibrium tethered 
shapes for each extension. However, even with the viscous regularization, the mesh can still be distorted by the 
dramatic deformations experienced at larger extensions. Therefore, re-meshing is performed at intervals of the 
extension. Fig. 4 shows snapshots from a typical simulation for reduced volume v = 0.9. 
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Figure 4: The tethering of a vesicle starting from a prolate of reduced volume v =■ 0.9. Left: limit surfaces of 
equilibrium shapes; right: control meshes. Number of nodes and elements: (a) 4202 nodes, 8400 elements; (b) & 
(c) 4682 nodes, 9360 elements; (d) 6202 nodes, 12400 elements; (e) & (f) 8682 nodes, 17360 elements. End to 
end distance: (a) 6.8/im, (b) 8.2/im, (c) 9.2/im, (d) 10.2/im, (e) 11.6/im, and (f) 12.8/^m 
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Force vs. end-to-end distance 



End-to-end distance (/im) 


6.8 


8.2 


9.2 


10.2 


11.6 


12.8 


Computed tether radius (^m) 


n/a 


n/a 


0.20 


0.165 


0.140 


0.105 


Computed force (pN) 





1.41 


1.76 


2.14 


2.72 


3.68 


Analytical estimated force (pN) 


n/a 


n/a 


1.88 


2.29 


2.69 


3.59 


Computed tension (pN//.jm) 


0.05 


0.43 


0.71 


0.98 


1.60 


2.93 


Analytical estimated tension (pN//im) 


n/a 


n/a 


0.75 


1.10 


1.53 


2.72 



Table 3: Computed and analytical estimated forces for each tethered shape. For the shape in Fig. 4(c), the 
computed tether radius is not uniform along the extension direction, and r = 0.20/im is an estimate. As noted in 
the text a bending modulus of JCc = ISA'^T is used for all calculations. 



Applied forces. The reaction forces conjugate to specified end displacements can also be calculated by simply 
adding up all internal forces of the fixed nodes (Eq. 13). The force vs. end-to-end distance results for the vesicles 
in Fig. 4 are shown in Table 3, with r the radius of the tethers (/-tm) and bending modulus ICc = l5Ki,T [37]. 
Although an exact analytical solution for the force-extension relation is not possible, a simple analytical estimate 
[39] is used to compare with the computed results from the simulation. The estimate assumes that the thin tube 
(tether) is pulled out from a sphere, and the sphere remain unchanged during the pulling (Fig. 5). The analytical 
estimated force and surface tension are given as [39]: 

F = 27r/Cc/r, (23) 

and the surface tension 

a = 0.5/Cc/r^ (24) 




Figure 5: Schematic of the tethered shape 

As Table 3 shows, for well developed tethered shapes (end-to-end distance 1 1.6 and 12.8 /im, vesicle (e) and 
(f) in Fig. 4), the computed results and analytical estimations are very close. It is a notable advantage that the 
present simulation framework is also capable of force-extension calculations for shapes that are not as simple 
as the schematic in Fig. 5. Although the present example is in fact axisymmetric, the algorithms are fully three- 
dimensional and can be applied to loadings and shapes lacking symmetry. 

3.3 Lipid Phase Separation 

Membranes formed from different lipids can separate into distinct domains (phases) according to their chemical 
properties, leading to the formation of buds [42, 43]. Baumgart et al. [44] found that their experiments are in 
good agreement with line tension theory [45-47], which treats domain interfaces as discrete with an interface 
energy proportional to their length. An alternative, smooth-interface approach, based on traditional Ginzburg- 
Landau (GL) theory [48, 49] can be used to also model phase separation [50-54]. One major drawback of the line 
tension model is that it requires the system to be pre-phase-separated into well-defined domains, preventing the 
consideration of composition dynamics. 

Here a GL model for a multi-component bilayer with two different lipids in equilibrium is formulated, assum- 
ing that the vesicle is composed of a mixture of two lipids denoted A and B. In general, these two lipid types 
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-0,5 03 1 1,5 
c 

Figure 6: Double-well normalized GL energy ip{c) ~ 16c^(c — 1)^ used to model phase segregation in a two- 
component lipid system. 

may have different constitutive properties, as modeled by separate constitutive parameters: {/c[^^ , /C^^'' , Cq'^^ } 
for lipid A, and {/c[f \ /c[f ^ , C^^^'' } for lipid B. 

Let the local concentrations of the two lipids be described by the concentration parameters c'"^-* , c*-^' G [0, 1] 
with c^^^ + c^^-* = 1. The local lipid concentration at point s = (s^, s^) can then be described by an order- 
parameter field c(s) = c^^\ which is referred to as the concentration field or phase field. The local constitutive 
properties of the membrane can then be modeled as functions of the phase field with convex combinations of the 
pure phase parameters; 









c)4^) 


(25a) 


JCg 








(25b) 


Co 








(25c) 



Thus rewriting the strain energy including explicit dependence of fields on surface position, 

E = |^/Cc(s)[2i/(s) - Co(s)]2 +/Cg(s)A'(s)| V^rf's (26) 

where explicit dependence of the mechanical properties on surface coordinates s has been noted as a reminder of 
the heterogeneity of the system. 

The mechanics of the membrane are then dependent on both the shape of the vesicle and the lipid compo- 
sition. Minimization of the total potential energy now yields two sets of Euler-Lagrange equations, one being 
the equilibrium equations related to variations in the shape 5x, and the other being a phase equilibrium equation 
related to variations in the concentration 5c. 

One further modification to the energy functional is needed to build into the model of the physics of phases 
separation [49]. 

I^E+ ( /:^E[i!{c)+e^\\7c\^]y/^(f-s. (27) 

JM 

Here the normalized GL energy ?/)(c) is a double- well potential such as 

ip{c) ^ 16c2(c- 1)2 

(see Fig. 6) which is minimized when the concentration c takes a value of either or 1, corresponding to local 
lipid concentration of either pure type A or pure type B. The parameter AE scales the height of the barrier 
between the two minima of tp{c), and controls the energy cost of a domain interface. The second addition to the 
energy describes short-range cooperativity between neighboring lipids. The parameter e is essentially a length 
scale which will determine the width of the region of transition between phases. As e decreases to zero, this 
region will Hmit to a curve where the concentration gradient can be non-zero. Inclusion of this penalty term in the 
energy will then produce the effect of a diffuse line tension in the transition between regions of pure phases. 

In Baumgart's experiment [42, 44], bending modulus ICc ~ 10~^^ J, line tension <t « 10~^'^N, and the 
radius of the vesicle R ~ lOfim. Two vesicles from [44] are simulated: one with reduced volume i' = 0.98, 
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phase B area fraction xb = (1/^) / = 0.89; the other v = 0.76, xb = (1/^) / = 0.56, starting from 
the original spherical shape with two separate domains (Fig. 2A and 2G in [44]). For the first one {i/ = 0.98), the 
simulation captured the small cap seen in the experiment (Fig. 7). 

But for 1/ = 0.76, the mesh around the interface is distorted (Fig. 8 (a)). Because the simulation starts 
from a sphere with roughly equilateral triangle elements, the shape change of vesicle causes the elements in 
the interface region to contract severely in the circumferential direction. This Element distortion needs to be 
suppressed because it can lead to inaccuracy and instability of the finite element simulation. 

In order to get a good mesh after the deformation, the elements near the interface need to contract in all 
directions so that they remain equilateral, resulting in a greater density of elements than other parts of the vesicle. 
To tackle such large deformations, remeshing strategies are often needed. The viscous regularization introduced 
in section 2.2 makes an r-adaptive remeshing possible, since reference configuration can be arbitrarily formulated 
to reposition the nodes of the mesh. Here, a slight modification of the dashpot regularization method is proposed 
with a reference updating strategy that drives elements toward equilateral shape. 

Given an element of the mesh at regularization iteration n — 1 with area An-i, r-adaptive regularization at 
step n is defined by placing springs on the three edges all of the same reference length 

i.e., the length of a side of an equilateral triangle of the same area An~i- Thus the regularization energy term for 
each triangle is written as 

,3 

En = ^Y.(ei,-in-l)\ (28) 

i=l 

where the are the lengths of the element edges. 

In principle this regularization energy could be applied to every element in a mesh. However, in practice these 
iterative updates are slow to converge to a fully relaxed state (with zero regularization energy), and depending on 
the regularization constant k the method can get stuck in a state with finite energy stored in the springs. Hence, 
this "equilateral" form of the dashpot regularization is only applied selectively to poorly shape elements, all the 
other elements with the standard viscous regularization. A shape-criteria 7 is then formulated to calculate different 
spring energy for different elements. 

Using this measure of shape quality, the regularization energy is defined for each triangle by 

£, ^ ^|ELl(C-4-l)^ 7 large, 
" UeLi(C-C-i)', 7small. 

In other words, if 7 is large (say, 7 > 1) for an element, it has poor shape and r-adaptive regularization is used 
on that element; if 7 is small enough, reference lengths are updated from the deformed lengths of the previous 
iteration as for the dashpot model described earlier. The addition of r-adaptive regularization has the effect of 
moving the nodes around on the membrane surface. In the present example of a phase-separated vesicle, this 
results in a finer mesh near the interface area than elsewhere (Fig. 8 (b)). For reduced volume i' = 0.76, the 
simulated result is shown in Fig. 7 to compare with the experimental result. 

4 Conclusions 

In this paper a framework for three-dimensional analysis of mechanics of lipid bilayer membranes is presented, 
based on the finite element method. Particular interest is focused on large deformation problems: tether formation 
(Sec. 3.2) and phase separation (Sec. 3.3). 

The primary difficulty faced in FE simulation of fluid membranes is the presense of mesh instabilities linked 
to the parameterization-independent nature of fluid surfaces. Curvature models of vesicle mechanics depend 
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(c) Simulation (d) Experiment (Fig. 2A from [44]) 



Figure 7: Comparison of simulation and experiment for two component lipid phase separation, (a) & (b) have 
reduced volume v = 0.98 and global concentration xb = 0.98. (c) & (d) have v = 0.76 and xb = 0.56. Phase 
A {c = 0) is colored blue; phase B (c = 1) is colored red. Scale bars are 5/^jm. Experimental images are taken 
from the work of Baumgart et al. [44]. 
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(a) standard dashpot regularization (b) r-adaptive regulaiization 



Figure 8: r-adaptive regularization helps elements perform well in the interface (ly = 0.76). 

only on current shape, and thus is not sensitive to in-plane (stretching and shearing) deformations of the surface 
FE mesh. Here a viscous regularization method is thus introduced to regularize tangential mesh deformations. 
In this method artificial reference configurations and corresponding in-plane energies are added to stabilize the 
tangential deformations; reference updates are designed so that artificial energy converges to zero in order to retain 
the physics of the original model. 

Regulariztion of tangential mesh deformations eliminates the need for local enforcement of membrane in- 
compressibility [17], providing a more convenient setting for augmented Lagrangian (AL) enforcement of global 
constraints on area and volume. The AL method can achieve higher accuracy with lower computational cost, 
compared to the penalty method. 

Large deformation problems can be very sensitive to mesh quality. Because of the physical meaninglessness 
of the reference configurations in the simulation, r-adaptive remeshing is easy to achieve in the context of viscous 
regularization, simply choosing a reference updating strategy which will reposition the nodes to get a better quality 
mesh. 

One promising direction for future work is to combine viscous regularization and r-adaptive remeshing with 
the dynamic triangulation approach [18], in which the edge of a pair of triangles swaps to form less distorted 
triangle elements instantly. This could be a powerful approach, speeding up the otherwise slow movements of 
nodes driven by the viscous regularization. Also, the success of r-adaptive regularization relies to some degree on 
the quality of the starting mesh. If there are too many badly shaped element and the shape criteria 7 tolerance is 
chosen to be too small, the mesh can sometimes lock with non-zero regularization energy, resulting in physically 
wrong shapes of vesicles. Dynamic triangulation could aleviate such locking. Lastly a r-adaptive regularized 
dynamic triangulation strategy could avoid the need for global remeshing in large deformation problems such as 
the tether simulations in Sec. 3.2. 

Although the problems simulated in this paper are all axisymmetric, the model is really designed for fully 
three-dimensional calculations, and can thus deal with arbitrary geometries and loads. For example future ap- 
plications such as mechanics of organelles like mitochondria [55] and endoplasmic reticulum (ER) [56], with 
incredibly complex shapes may provide exciting opportunities for future study with these methods. 
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